function [Tipw,std_T_ipw,std_b1ipw,std_b2ipw,loweripw,upperipw,...
          cri1ipwmed,cri2ipwmed,cri1ipw,cri2ipw,QTEModelSelIdx,ATECVSelectIdx] = ...
          CVSieveBasisIPWType2(Y,X,A,n,taulist,theta_true,B)
% This function obtains the CV IPW statistics 
% Inputs:
%       Y:          the targeted variables, 2n*1
%       X:          the associated basic covariats, 2n*K
%       A:          the contrl or treatment status, 2n*1
%       n:          the number of pairs
%       taulist:    the quantiles aim to examine, L*1
%       theta_true: the hypothesis test value 
%       B:          the number of Boostrap samples
% Outputs: 
% Note: Q(tau), is the tau-quantiles of Boostrap samples
%       Tipw:            the values of T-statistic for the hypothesis
%       std_T_ipw:       the standard errors of T-statistic 
%       std_b1ipw:       the standard error based Q(0.75) and Q(0.25)
%       std_b2ipw:       the standard error The standard error of CV IPW QTE based Q(0.975) and Q(0.025)
%       loweripw:        the Q(0.025) of Boostrap samples
%       upperipw:        the Q(0.975) of Boostrap samples
%       cri1ipwmed:      the critical values (C_α in Section 7) in the the 
%                        Boostrap methods given median Q(0.5) and std_b1breve
%       cri2ipwmed:      the critical values (C_α in Section 7) in the the 
%                        Boostrap methods given median Q(0.5) and std_b2breve
%       cri1ipw:         the critical values (C_α in Section 7) in the the 
%                        Boostrap methods given median (Q(0.975) + Q(0.025))/2 and std_b1breve
%       cri2ipw:         the critical values (C_α in Section 7) in the the 
%                        Boostrap methods given median (Q(0.975) + Q(0.025))/2 and std_b2breve
% Date: 2020-04-20

    % 1. The pre-run process of observations
    n1     = sum(A);
    n0     = 2*n - n1;

    Y1bar  = sum(Y.*A)/n1;
    Y0bar  = sum(Y.*(1-A))/n0;

    L      = length(taulist);
    % 2. The storages of the Boostrap samples for the boostrap
    q1ipw  = zeros(B,L);
    q0ipw  = zeros(B,L);
    qipw   = zeros(B,L + 3);


    % 3. The generating of Boostrap innovations
    N      = 2*n;
    xi     = -log(1-rand(B,N)); % weights for IPW

    % 4. The computation of the sieve basis weights when dim(X) = 1 or 2.
    [Ahatmat, thetaipw, QTEModelSelIdx, ATECVSelectIdx,xi] = CVScoreCalType2(Y,X,A,taulist,xi,B);

    % 5. Generating Boostrap samples for the quantiles given in the
    % taulist.
    for t = 1:length(taulist)
        tau           = taulist(t);
        % The sieve basis weighted IPW estimator
        q1ipw(:,t)    = myqr1Matrix(Y,A,Ahatmat(:,:,t),tau*N,xi);
        q0ipw(:,t)    = myqr0Matrix(Y,A,Ahatmat(:,:,t),tau*N,xi);
        qipw(:,t)     = q1ipw(:,t)-q0ipw(:,t); 
    end

    % The sieve basis weighted IPW ATE and IPW QTE 
    std_T_ipw   = ((prctile(thetaipw,97.5) - prctile(thetaipw,2.5)))/(norminv(0.975) - norminv(0.025));
    Tipw        = (Y1bar - Y0bar - theta_true)/std_T_ipw;
    % Compute the case Q(tau1) - Q(tau2)
    qipw(:,L+1) = qipw(:,L-1) - qipw(:,1); % 75% quantile - 25% quantile
    qipw(:,L+2) = qipw(:,L) - qipw(:,1); % 50% quantile - 25% quantile
    qipw(:,L+3) = qipw(:,L-1) - qipw(:,L);  % 75% quantile - 50% quantile

    % 5. The computation of the two types of standard errors:
    %       - sig1 is based on Q(0.75) and Q(0.25)
    %       - sig2 is based on Q(0.975) and Q(0.025)
    % and the Boostrap samples quantiles Q(0.025) and Q(0.975)    
    std_b1ipw   = ((prctile(qipw,75) - prctile(qipw,25)))/(norminv(0.75) - norminv(0.25));
    std_b2ipw   = ((prctile(qipw,97.5) - prctile(qipw,2.5)))/(norminv(0.975) - norminv(0.025));    
    loweripw    = prctile(qipw,2.5);
    upperipw    = prctile(qipw,97.5);

    % 6. The computation of critical values (C_α in Section 7) in the 
    % the Boostrap methods given median Q(0.5) and (Q(0.975) + Q(0.025))/2
    % and based on the two standard errors above.

    % The median Q(0.5) case
    cri1ipwmed  = prctile(max( ...
        abs(qipw(:,1:L)-repmat(median(qipw(:,1:L)),B,1))./repmat(squeeze(std_b1ipw(1:L)),B,1),...
        [],2),95);
    cri2ipwmed  = prctile(max( ...
        abs(qipw(:,1:L)-repmat(median(qipw(:,1:L)),B,1))./repmat(squeeze(std_b2ipw(1:L)),B,1),...
        [],2),95);    

    % The (Q(0.975) + Q(0.025))/2 case
    cri1ipw     = prctile(max( ...
        abs(qipw(:,1:L)-repmat((prctile(qipw(:,1:L),97.5) + prctile(qipw(:,1:L),2.5))/2,B,1))./repmat(squeeze(std_b1ipw(1:L)),B,1),...
        [],2),95);
    cri2ipw     = prctile(max( ...
        abs(qipw(:,1:L)-repmat((prctile(qipw(:,1:L),97.5) + prctile(qipw(:,1:L),2.5))/2,B,1))./repmat(squeeze(std_b2ipw(1:L)),B,1),...
        [],2),95);
    
end

% Quantile Treatment Effects and Bootstrap Inference under Covariate-Adaptive Randomization 
% Yichong Zhang, Xin Zheng
% 2020/02/14

% [Generate q0 in IPW estimation] 

% Note: It is different from theirs. Here W is of size B*N, B is the number
% of Boostrap samples and N is the number observations

function q0 = myqr0Matrix(Y,T,P,k,W)
    [B,n]       = size(W);
    [eY,idx]    = sort(Y.*(1-T)+(10000+Y).*T);
    L           = length(k);
    if L ~= 1
        error('myqr0Matrix error');
    end
    
    
    if size(P,2) ~=1
        k1          = sum(repmat(1 - T',B,1).*W./(1-P + eps),2)*k/n;

        Wtranspose = W';
        idxMatrix    = repmat(idx,1,B) + repmat(0:B-1,n,1) * n;
        part01       = Wtranspose(idxMatrix);
        Ptranspose = P';
        e          = part01./(1-Ptranspose(idxMatrix)+eps);

        cumsumMatrixe01   = cumsum(e);
        cumsumMatrixe01(1,:) = 0;
        cumsumMatrixe02   = nan(size(cumsumMatrixe01));
        cumsumMatrixe02(2:end,:) = cumsumMatrixe01(1:end-1,:);

        IterStopIdx       = floor(k1) + 1;

        TmpLessIdx      = (cumsumMatrixe01 >= repmat(k1',n,1));
        TmpLargeIdx     = (cumsumMatrixe02 < repmat(k1',n,1));
        FinalIdx        = TmpLargeIdx.*TmpLessIdx;
        FinalIdx(cumsum(FinalIdx) > 1)    = 0;
        [I,J]           = find(FinalIdx == 1);
        FinalIdx(:,J(I > IterStopIdx(J))) = 0;


        EmatrixFinal    = repmat(eY,1,B).*FinalIdx;
        EmatrixFinal(EmatrixFinal == 0) = nan;
        q0              = max(EmatrixFinal);
        q0(isnan(q0))   = 0;     
        
    else
        k1          = sum(repmat(1 - T',B,1).*W./repmat(1 - P'+eps',B,1),2)*k/n;

        Wtranspose  = W';
        part01      = Wtranspose(repmat(idx,1,B) + repmat(0:B-1,n,1) * n);
        e           = part01./(P(repmat(idx,1,B))+eps);
        cumsumMatrixe01   = cumsum(e);
        cumsumMatrixe01(1,:) = 0;
        cumsumMatrixe02   = nan(size(cumsumMatrixe01));
        cumsumMatrixe02(2:end,:) = cumsumMatrixe01(1:end-1,:);

        IterStopIdx       = floor(k1) + 1;

        TmpLessIdx      = (cumsumMatrixe01 >= repmat(k1',n,1));
        TmpLargeIdx     = (cumsumMatrixe02 < repmat(k1',n,1));
        FinalIdx        = TmpLargeIdx.*TmpLessIdx;
        [I,J]           = find(FinalIdx == 1);
        FinalIdx(:,J(I > IterStopIdx(J))) = 0;

        EmatrixFinal    = repmat(eY,1,B).*FinalIdx;
        EmatrixFinal(EmatrixFinal == 0) = nan;
        q0              = max(EmatrixFinal);
        q0(isnan(q0))   = 0;
    end
end


% Quantile Treatment Effects and Bootstrap Inference under Covariate-Adaptive Randomization 
% Yichong Zhang, Xin Zheng
% 20-04-13
% [Generate q1 in IPW estimation]
% Note: It is different from theirs. Here W is of size B*N, B is the number
% of Boostrap samples and N is the number observations

function q1 = myqr1Matrix(Y,T,P,k,W)
% Here k must be 1-dimensional
    [B,n]       = size(W);
    [eY,idx]    = sort(Y.*T+(10000+Y).*(1-T));
    L           = length(k);
    if L ~= 1
        error('myqr1Matrix error');
    end
    
    
    if size(P,2) ~= 1
        k1    = sum(repmat(T',B,1).*W./(P+eps),2)*k/n;
        Wtranspose = W';
        idxMatrix    = repmat(idx,1,B) + repmat(0:B-1,n,1) * n;
        part01       = Wtranspose(idxMatrix);
        Ptranspose = P';
        e     = part01./(Ptranspose(idxMatrix)+eps);
        %e(:,myqrb) - e(:,2)

        cumsumMatrixe01   = cumsum(e);
        cumsumMatrixe01(1,:) = 0;
        cumsumMatrixe02   = nan(size(cumsumMatrixe01));
        cumsumMatrixe02(2:end,:) = cumsumMatrixe01(1:end-1,:);

        IterStopIdx       = floor(k1) + 1;

        TmpLessIdx      = (cumsumMatrixe01 >= repmat(k1',n,1));
        TmpLargeIdx     = (cumsumMatrixe02 < repmat(k1',n,1));
        FinalIdx        = TmpLargeIdx.*TmpLessIdx;
        FinalIdx(cumsum(FinalIdx) > 1)    = 0;
        [I,J]           = find(FinalIdx == 1);
        FinalIdx(:,J(I > IterStopIdx(J))) = 0;

        EmatrixFinal    = repmat(eY,1,B).*FinalIdx;
        EmatrixFinal(EmatrixFinal == 0) = nan;
        q1              = max(EmatrixFinal);
        q1(isnan(q1))   = 0;
    else
        k1          = sum(repmat(T',B,1).*W./repmat(P'+eps',B,1),2)*k/n;

        Wtranspose  = W';
        part01      = Wtranspose(repmat(idx,1,B) + repmat(0:B-1,n,1) * n);
        e           = part01./(P(repmat(idx,1,B))+eps);
        cumsumMatrixe01   = cumsum(e);
        cumsumMatrixe01(1,:) = 0;
        cumsumMatrixe02   = nan(size(cumsumMatrixe01));
        cumsumMatrixe02(2:end,:) = cumsumMatrixe01(1:end-1,:);

        IterStopIdx       = floor(k1) + 1;

        TmpLessIdx      = (cumsumMatrixe01 >= repmat(k1',n,1));
        TmpLargeIdx     = (cumsumMatrixe02 < repmat(k1',n,1));
        FinalIdx        = TmpLargeIdx.*TmpLessIdx;
        FinalIdx(cumsum(FinalIdx) > 1)    = 0;
        [I,J]           = find(FinalIdx == 1);
        FinalIdx(:,J(I > IterStopIdx(J))) = 0;
        

        EmatrixFinal    = repmat(eY,1,B).*FinalIdx;
        EmatrixFinal(EmatrixFinal == 0) = nan;
        q1              = max(EmatrixFinal);
        q1(isnan(q1))   = 0;
        
    end
end







